{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "75e2c1fe-8a06-4f0b-8fd6-d1b6d10ce82f",
   "metadata": {},
   "source": [
    "# Quantum Phase Estimation for Solving Matrix Eigenvalues\n",
    "\n",
    "Quantum Phase Estimation (QPE) is a key algorithm in quantum computing, allowing you to estimate the phase (or eigenvalue) relating to a Hermitian matrix. The algorithm is designed such that given the inputs of a matrix $M$ and an eigenvalue ${|\\psi\\rangle}$, the output obtained is $\\theta$, where\n",
    "\n",
    "$ U{|\\psi\\rangle} = e^{2\\pi i\\theta}{|\\psi\\rangle} , U = e^{2\\pi iM} $.\n",
    "\n",
    "By measuring the accumulated phase, the QPE algorithm calculates the eigenvalues relating to the chosen input vector. To read more about the QPE algorithm and its method for achieving the phase, refer to [[1](#NC)].\n",
    "\n",
    "Generally speaking, when the eigenvectors of the matrix are not known in advance yet the eigenvalues are sought, you can choose a random vector ${|v\\rangle}$ for the algorithm\u2019s initial state. Some eigenvalues will be found as the vector can be described in the matrix's basis, defined by the set of eigenvalues of $M$: {$\\psi_i$}. Generally, any vector can be written as a superposition of any basis set, thus\n",
    "\n",
    "${|v\\rangle} = \\sum_i a_i{|\\psi_i\\rangle}$\n",
    "\n",
    "and\n",
    "\n",
    "$U{|v\\rangle} = \\sum_i a_i e^{2\\pi i\\theta_i}{|\\psi_i\\rangle}$.\n",
    "\n",
    "Using execution with enough shots, you can obtain this set of $\\theta_i$; i.e., a subset of the matrix's eigenvalues.\n",
    "\n",
    "**This tutorial presents a generic usage of the QPE algorithm:**\n",
    "\n",
    "1. Define a matrix.\n",
    "\n",
    "2. Initialize a state either with its eigenstate or with a random vector.\n",
    "\n",
    "3. Choose a resolution for the solution.\n",
    "\n",
    "4. Find the related eigenvalues using QPE and analyze the results."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cd71faa-1ef3-4d67-8c9a-d66e0a895b9a",
   "metadata": {},
   "source": [
    "## Prerequisites\n",
    "\n",
    "This tutorial uses external libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "11eaf8f8-c06e-443e-ac74-5aa364f5198c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:09.334719Z",
     "iopub.status.busy": "2024-05-07T14:54:09.334236Z",
     "iopub.status.idle": "2024-05-07T14:54:09.467189Z",
     "shell.execute_reply": "2024-05-07T14:54:09.458037Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import itertools  # noqa\n",
    "import math\n",
    "from itertools import product\n",
    "from typing import List, cast\n",
    "\n",
    "import numpy as np\n",
    "from numpy import kron, linalg as LA"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac537458-568d-4718-9020-408cc7232641",
   "metadata": {},
   "source": [
    "## 1. Setting a Specific Example"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a6f04f8-f4f7-46dc-a8e3-05559d8906f3",
   "metadata": {},
   "source": [
    "### 1.1. Set the Matrix\n",
    "\n",
    "Define the matrix to submit. This can be any Hermitian matrix with size $2^n$ by $2^n$ with $n$ a positive integer. Throughout the code this matrix is given in the variable `M`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3e9c8c83-58ed-4f99-89ac-9cea3f1b4a5e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:09.471500Z",
     "iopub.status.busy": "2024-05-07T14:54:09.471139Z",
     "iopub.status.idle": "2024-05-07T14:54:09.478144Z",
     "shell.execute_reply": "2024-05-07T14:54:09.477224Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "M = np.array([[0, 3, 4, 0], [-0.8, 3, 0, 0], [1, 0, -0.5, 5], [0, 0, 0, -0.75]])\n",
    "\n",
    "M = np.array(\n",
    "    [\n",
    "        [0.38891555, 0.23315811, 0.21499372, 0.06119557],\n",
    "        [0.23315811, 0.44435328, 0.25197881, -0.13087919],\n",
    "        [0.21499372, 0.25197881, 0.44116509, -0.01961855],\n",
    "        [0.06119557, -0.13087919, -0.01961855, 0.32556608],\n",
    "    ]\n",
    ")\n",
    "\n",
    "M_t = M.transpose()\n",
    "\n",
    "M = (M + M_t) / 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a6edbc6-1100-4824-a15e-c985728f544b",
   "metadata": {},
   "source": [
    "### 1.2. Set the Initial Vector \n",
    "\n",
    "Choose the vector that will be defined later as the initial condition for the run. There are two options:\n",
    "1. Define your own initial vector in the variable `int_vec`, while setting the parameter `eigen_vec` as `False`.\n",
    "2. Set `eigen_vec` to `True`, then you can choose the index `ev` of the eigenvalue that will be set as the initial state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "736b0585-d019-4b22-84af-1e7f36a509d5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:09.481875Z",
     "iopub.status.busy": "2024-05-07T14:54:09.481400Z",
     "iopub.status.idle": "2024-05-07T14:54:09.487794Z",
     "shell.execute_reply": "2024-05-07T14:54:09.487109Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Your initial state is [0.80216327 0.54823687 0.6653769  0.62012979]\n"
     ]
    }
   ],
   "source": [
    "eigen_vec = False\n",
    "ev = 1\n",
    "\n",
    "if eigen_vec:\n",
    "    w, v = LA.eig(M)\n",
    "    print(\"the eigenvalues are\", w)\n",
    "    print(\"the eigenvectors are\", v, sep=\"\\n\")\n",
    "    int_vec = v[:, ev]\n",
    "else:\n",
    "    int_vec = np.random.rand(np.shape(M)[0])\n",
    "\n",
    "print(\"Your initial state is\", int_vec)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7b4d3fb-f5cc-45d0-b02c-d74d74677de3",
   "metadata": {},
   "source": [
    "## 2. Constructing Auxiliary Functions\n",
    "\n",
    "Defining some auxiliary functions is essential for designing the QPE in a modular fashion."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7785e9ca-3532-485c-80f8-e405a4a419ab",
   "metadata": {},
   "source": [
    "### 2.1 Classical Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6876d85-70ca-4be7-bf2e-6a78235a77e5",
   "metadata": {},
   "source": [
    "#### 2.1.1 Pauli Decomposition\n",
    "\n",
    "To translate the matrix into quantum circuit language, write the matrix in the form of a list of strings of composite Pauli operators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3296cfde-4155-4d76-9ba6-16aa24677a61",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:09.492717Z",
     "iopub.status.busy": "2024-05-07T14:54:09.492293Z",
     "iopub.status.idle": "2024-05-07T14:54:09.504049Z",
     "shell.execute_reply": "2024-05-07T14:54:09.503371Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "Paulidict = {\n",
    "    \"I\": np.array([[1, 0], [0, 1]], dtype=np.complex128),\n",
    "    \"Z\": np.array([[1, 0], [0, -1]], dtype=np.complex128),\n",
    "    \"X\": np.array([[0, 1], [1, 0]], dtype=np.complex128),\n",
    "    \"Y\": np.array([[0, -1j], [1j, 0]], dtype=np.complex128),\n",
    "}\n",
    "\n",
    "\n",
    "# generate all combinations of Pauli strings of size n\n",
    "def generate_all_pauli_strings(seq, n):\n",
    "    for s in product(seq, repeat=n):\n",
    "        yield \"\".join(s)\n",
    "\n",
    "\n",
    "# convert a Paulistring of size n to 2**n X 2**n matrix\n",
    "def pauli_string_2mat(seq):\n",
    "    myPmat = Paulidict[seq[0]]\n",
    "    for p in seq[1:]:\n",
    "        myPmat = kron(myPmat, Paulidict[p])\n",
    "    return myPmat\n",
    "\n",
    "\n",
    "# Hilbert-Schmidt-Product of two matrices M1, M2\n",
    "def hilbert_schmidt(M1, M2):\n",
    "    return (np.dot(M1.conjugate().transpose(), M2)).trace()\n",
    "\n",
    "\n",
    "# naive decomposition, running over all HS products for all Pauli strings\n",
    "def lcu_naive(H):\n",
    "    assert H.shape[0] == H.shape[1], \"matrix is not square\"\n",
    "    assert H.shape[0] != 0, \"matrix is of size 0\"\n",
    "    assert H.shape[0] & (H.shape[0] - 1) == 0, \"matrix size is not 2**n\"\n",
    "\n",
    "    n = int(np.log2(H.shape[0]))\n",
    "    myPauliList = list(generate_all_pauli_strings(\"IZXY\", n))\n",
    "\n",
    "    mylist = []\n",
    "\n",
    "    for pstr in myPauliList:\n",
    "        co = (1 / 2**n) * hilbert_schmidt(pauli_string_2mat(pstr), H)\n",
    "        if co != 0:\n",
    "            mylist = mylist + [(pstr, co)]\n",
    "\n",
    "    return mylist"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9303ffc8-9579-49e7-a498-7899d2bf1624",
   "metadata": {},
   "source": [
    "#### 2.1.2 Parser to the `PauliTerm` struct\n",
    "\n",
    "A Hamiltonian is defined by a list of PauliTerm structs. Define a function that transforms between the Pauli list representation to this struct presentation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "90918038-cd4c-47c8-b8e5-b026315e5874",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:09.508811Z",
     "iopub.status.busy": "2024-05-07T14:54:09.507613Z",
     "iopub.status.idle": "2024-05-07T14:54:12.125675Z",
     "shell.execute_reply": "2024-05-07T14:54:12.122731Z"
    }
   },
   "outputs": [],
   "source": [
    "from classiq import Pauli, PauliTerm\n",
    "\n",
    "my_list = {\"I\": Pauli.I, \"X\": Pauli.X, \"Y\": Pauli.Y, \"Z\": Pauli.Z}\n",
    "\n",
    "\n",
    "def pauli_str_to_enums(pauli):\n",
    "    return [my_list[s] for s in pauli]\n",
    "\n",
    "\n",
    "def pauli_operator_to_hamiltonian(pauli_list):\n",
    "    return [\n",
    "        PauliTerm(\n",
    "            pauli=pauli_str_to_enums(pauli), coefficient=cast(complex, coeff).real\n",
    "        )\n",
    "        for pauli, coeff in pauli_list\n",
    "    ]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dc21993-ef3e-48af-b831-ac8fd101067f",
   "metadata": {},
   "source": [
    "#### 2.1.3 Matrix Rescaling"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f1896e5-fda6-48e3-b12c-259670c5b8d8",
   "metadata": {},
   "source": [
    "As QPE obtains a phase in the form $e^{2\\pi i\\theta}$, there is meaning only for $\\theta \\in [0,1)$. However, the matrix M can have any eigenvalue. To fix this discrepancy, the values of the matrix stretch to be rescaled. If\n",
    "$\\theta \\in [\\lambda_{min}, \\lambda_{max}]$ you can use a normalization function to map those values into $[0, 1-1/{2^m}]$, where $m$ \n",
    "is the size of the QPE register.\n",
    "\n",
    "Perform the normalization procedure by:\n",
    "\n",
    "a. Defining the function `normalization_params()` that finds a rough estimation for the eigenvalue with the largest absolute value by adding together all the Pauli coefficients and multiplying by the matrix's dimensions. This yields a value $\\lambda$ (which is referred to in the code as `normalization_coeff`) and now you can assume that the domain is $\\theta \\in [-\\lambda, \\lambda]$.\n",
    "In general, you can build a more accurate assessment that decreases the span of solutions and thus achieves a better resolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c4a4efd0-06e2-461f-8a37-14aeaa4cf1fb",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:12.132642Z",
     "iopub.status.busy": "2024-05-07T14:54:12.129665Z",
     "iopub.status.idle": "2024-05-07T14:54:12.138315Z",
     "shell.execute_reply": "2024-05-07T14:54:12.135759Z"
    }
   },
   "outputs": [],
   "source": [
    "def normalization_params(pauli_list):\n",
    "    return len(pauli_list[0][0]) * sum(\n",
    "        [abs(pauli_list[k][1]) for k in range(len(pauli_list))]\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d6c378f-ae23-4baf-999d-e9943a7e49b3",
   "metadata": {},
   "source": [
    "b. Defining the function `normalize_hamiltonian` that shifts the matrix by adding $\\lambda*I^n$ to the Pauli list. (The evaluated span is thus $\\theta\\in[0, 2*\\lambda]$) and normalizes it by multiplying all the Pauli coefficients by $(1-1/2^n)/(2*\\lambda)$ (the evaluated span is then $\\theta\\in [0, 1-1/2^n]$, as required.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "15ca935c-04a3-40c6-9846-f1b822bc9832",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:12.141568Z",
     "iopub.status.busy": "2024-05-07T14:54:12.141308Z",
     "iopub.status.idle": "2024-05-07T14:54:12.154818Z",
     "shell.execute_reply": "2024-05-07T14:54:12.154042Z"
    }
   },
   "outputs": [],
   "source": [
    "def normalize_hamiltonian(pauli_list, normalization_coeff, k):\n",
    "    list_size = len(pauli_list)\n",
    "    num_qubits = len(pauli_list[0][0])\n",
    "    normalization = (1 - 1 / (2**k)) / (2 * normalization_coeff)\n",
    "    normalized_list = [\n",
    "        (pauli_list[k][0], pauli_list[k][1] * normalization) for k in range(list_size)\n",
    "    ]\n",
    "    if \"I\" * num_qubits in [pauli_list[k][0] for k in range(list_size)]:\n",
    "        id_index = [y[0] for y in pauli_list].index(\"I\" * num_qubits)\n",
    "        normalized_list[id_index] = (\n",
    "            \"I\" * num_qubits,\n",
    "            (pauli_list[id_index][1] + normalization_coeff) * normalization,\n",
    "        )\n",
    "    else:\n",
    "        normalized_list.append((\"I\" * num_qubits, normalization_coeff * normalization))\n",
    "\n",
    "    return normalized_list"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c2cf50a-ef8a-41d2-8684-61c242734ed4",
   "metadata": {},
   "source": [
    "#### 2.1.3 QPE Precision Estimator"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2e8d8f4-4ff7-4444-a9dc-1abb7f36e653",
   "metadata": {},
   "source": [
    "For QPE algorithms, the precision is set by phase register size $m$, such that the resolution is $1/{2^m}$. If the matrix needs to be normalized, the resolution will be distorted. In the case of normalization, the span of results for the QPE stretches between the lowest and highest possible phase, thus the resolution is mapped to $normalization-coefficient/{2^m} ~\\sim 1/{((\\lambda_{max}-\\lambda_{min})*2^m)}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "09d5727b-25f3-423b-bfe4-03443d596877",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:12.157983Z",
     "iopub.status.busy": "2024-05-07T14:54:12.157724Z",
     "iopub.status.idle": "2024-05-07T14:54:12.163621Z",
     "shell.execute_reply": "2024-05-07T14:54:12.162868Z"
    }
   },
   "outputs": [],
   "source": [
    "def get_qpe_precision(pauli_list, desired_resolution):\n",
    "    nqpe = math.log2(2 * normalization_params(pauli_list) / desired_resolution)\n",
    "    return math.ceil(nqpe)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcb881a7-61ef-4e1d-b971-3fda5b1b0a78",
   "metadata": {},
   "source": [
    "## 2.2 Quantum Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db33c350-dfc6-42de-809a-9825fc055d7c",
   "metadata": {},
   "source": [
    "### 2.2.1 A Flexible QPE\n",
    "\n",
    "Define a flexible QPE function, which allows you to prescribe the \"telescopic\" expansion of the powered unitary via the `unitary_with_power` \"QCallable\" See the High Level Modeling Flexible QPE tutorial for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ba7a59fa-689c-4c95-9a4e-24064a46f98a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:12.166891Z",
     "iopub.status.busy": "2024-05-07T14:54:12.166652Z",
     "iopub.status.idle": "2024-05-07T14:54:12.175428Z",
     "shell.execute_reply": "2024-05-07T14:54:12.174663Z"
    }
   },
   "outputs": [],
   "source": [
    "from classiq import (\n",
    "    CInt,\n",
    "    H,\n",
    "    Output,\n",
    "    QArray,\n",
    "    QBit,\n",
    "    QCallable,\n",
    "    apply_to_all,\n",
    "    bind,\n",
    "    control,\n",
    "    invert,\n",
    "    qft,\n",
    "    qfunc,\n",
    "    repeat,\n",
    ")\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def my_qpe_flexible(\n",
    "    unitary_with_power: QCallable[CInt, QArray[QBit]],\n",
    "    precision: CInt,\n",
    "    state: QArray[QBit],\n",
    "    phase: Output[QArray[QBit, \"precision\"]],\n",
    ") -> None:\n",
    "    allocate(precision, phase)\n",
    "    apply_to_all(H, phase)\n",
    "\n",
    "    repeat(\n",
    "        count=precision,\n",
    "        iteration=lambda index: control(\n",
    "            ctrl=phase[index],\n",
    "            operand=lambda: unitary_with_power(2**index, state),\n",
    "        ),\n",
    "    )\n",
    "\n",
    "    invert(lambda: qft(phase))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad8f766e-91ed-44c5-9c4e-3f909358cdb1",
   "metadata": {},
   "source": [
    "### 2.2.2 A First Order Suzuki Trotter with power-logic\n",
    "\n",
    "Wrap the Trotter-Suzuki function of order 1 with a \"power-logic\" for the repetition as a function of its power."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "920ac18c-591a-4ec2-b178-451f173dd5da",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:12.178856Z",
     "iopub.status.busy": "2024-05-07T14:54:12.178501Z",
     "iopub.status.idle": "2024-05-07T14:54:12.185769Z",
     "shell.execute_reply": "2024-05-07T14:54:12.185021Z"
    }
   },
   "outputs": [],
   "source": [
    "from classiq import CArray, CReal, suzuki_trotter\n",
    "from classiq.qmod.symbolic import ceiling, log\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def suzuki_trotter1_with_power_logic(\n",
    "    hamiltonian: CArray[PauliTerm],\n",
    "    pw: CInt,\n",
    "    r0: CInt,\n",
    "    reps_scaling_factor: CReal,\n",
    "    evolution_coefficient: CReal,\n",
    "    target: QArray[QBit],\n",
    ") -> None:\n",
    "    suzuki_trotter(\n",
    "        hamiltonian,\n",
    "        evolution_coefficient=evolution_coefficient * pw,\n",
    "        order=1,\n",
    "        repetitions=r0 * ceiling(reps_scaling_factor ** (log(pw, 2))),\n",
    "        qbv=target,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13398c9e-a74d-4819-9205-158f00498e17",
   "metadata": {},
   "source": [
    "### 2.2.3 A Unitary with power-logic\n",
    "\n",
    "As an alternative to the Trotter-Suzuki formula, you can work with an exact unitary decomposition. In this case, the power-logic is naive:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "1496c3fb-d68d-4466-bfe1-828903f44218",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:12.188744Z",
     "iopub.status.busy": "2024-05-07T14:54:12.188485Z",
     "iopub.status.idle": "2024-05-07T14:54:12.199694Z",
     "shell.execute_reply": "2024-05-07T14:54:12.198097Z"
    }
   },
   "outputs": [],
   "source": [
    "from classiq import power, unitary\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def unitary_with_power_logic(\n",
    "    pw: CInt, matrix: CArray[CArray[CReal]], target: QArray[QBit]\n",
    ") -> None:\n",
    "    power(pw, lambda: unitary(elements=matrix, target=target))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffb62b82-d3ec-4e63-a558-b57af0c2e0b4",
   "metadata": {},
   "source": [
    "## 3. Preparing the Matrix for QPE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "9c95c12b-5783-456f-92ca-937df78b3259",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:12.202931Z",
     "iopub.status.busy": "2024-05-07T14:54:12.202668Z",
     "iopub.status.idle": "2024-05-07T14:54:12.210655Z",
     "shell.execute_reply": "2024-05-07T14:54:12.210012Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('II', (0.39999999999999997+0j)), ('IZ', (0.015040319999999996+0j)), ('IX', (0.10676978000000001+0j)), ('ZI', (0.016634415000000013+0j)), ('ZZ', (-0.042759185000000005+0j)), ('ZX', (0.12638833+0j)), ('XI', (0.042057264999999996+0j)), ('XZ', (0.172936455+0j)), ('XX', (0.15658719+0j)), ('YY', (0.09539162+0j))]\n"
     ]
    }
   ],
   "source": [
    "pauli_ops = lcu_naive(M)\n",
    "print(pauli_ops)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "1cc56c7a-8f2e-4f89-bc0c-7d768b0388d4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:12.214071Z",
     "iopub.status.busy": "2024-05-07T14:54:12.213616Z",
     "iopub.status.idle": "2024-05-07T14:54:12.219726Z",
     "shell.execute_reply": "2024-05-07T14:54:12.219064Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of qubits:  2\n"
     ]
    }
   ],
   "source": [
    "N = len(pauli_ops[0][0])\n",
    "print(\"number of qubits: \", N)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c69c18c8-566e-46e9-8818-b3baf69ff9eb",
   "metadata": {},
   "source": [
    "### 3.1 Choose the Algorithm's Precision\n",
    "\n",
    "Choose the precision using the `n_qpe` parameter or set your desired resolution. If you choose the resolution and set the parameter `get_recommended_qpe_size` parameter to True, the number of qubits is calculated for you accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "21c12d6d-399b-4ffe-9150-82adff52dd61",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:12.223259Z",
     "iopub.status.busy": "2024-05-07T14:54:12.222786Z",
     "iopub.status.idle": "2024-05-07T14:54:12.229367Z",
     "shell.execute_reply": "2024-05-07T14:54:12.228669Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of qubits for QPE is 8\n"
     ]
    }
   ],
   "source": [
    "n_qpe = 8\n",
    "\n",
    "# recommended QPE_SIZE:\n",
    "get_recommended_qpe_size = False\n",
    "\n",
    "desired_resolution = 0.02\n",
    "\n",
    "\n",
    "if get_recommended_qpe_size:\n",
    "    n_qpe = get_qpe_precision(pauli_ops, desired_resolution)\n",
    "\n",
    "print(\"number of qubits for QPE is\", n_qpe)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa3d1079-c12b-48db-8935-1ea059fdbce9",
   "metadata": {},
   "source": [
    "### 3.2 Normalize the Matrix\n",
    "\n",
    "Transform the matrix to ensure its eigenvalues are between $0$ to $1-(1/2^m)$. The QPE procedure is performed on the new normalized matrix. After the phases are obtained, gather the original phases of the pre-normalized matrix by performing opposite steps to this normalization procedure.\n",
    "\n",
    "* If the matrix eigenvalues are naturally between the values $0$ to $1-(1/2^n)$, you may not want to normalize them as that  may enlarge the span, thus lowering the resolution of the algorithm. In this case, skip those lines or change the value of `normalize` to False."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "5e986c5a-1aae-4d97-be9c-693d7d8c2ae0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:12.232828Z",
     "iopub.status.busy": "2024-05-07T14:54:12.232361Z",
     "iopub.status.idle": "2024-05-07T14:54:12.240203Z",
     "shell.execute_reply": "2024-05-07T14:54:12.239522Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('II', (0.582852238955473+0j)), ('IZ', (0.003188749529016948+0j)), ('IX', (0.02263662513086446+0j)), ('ZI', (0.0035267190456534513+0j)), ('ZZ', (-0.009065520615911005+0j)), ('ZX', (0.026796020813436065+0j)), ('XI', (0.00891670416324194+0j)), ('XZ', (0.036664847518610696+0j)), ('XX', (0.033198584096787005+0j)), ('YY', (0.020224302631005442+0j))]\n"
     ]
    }
   ],
   "source": [
    "# normalizing the operator\n",
    "## create a matrix such that its normalized version has eigenvalues of [0,1/2^k] where k is the resolution of the QPE\n",
    "normalize = True\n",
    "if normalize:\n",
    "    normalization_coeff = normalization_params(pauli_ops)\n",
    "    new_pauli_list = normalize_hamiltonian(pauli_ops, normalization_coeff, n_qpe)\n",
    "    print(new_pauli_list)\n",
    "\n",
    "    size = math.sqrt(M.size)\n",
    "    I = np.eye(int(size))\n",
    "\n",
    "    Mnew = (\n",
    "        (M + normalization_coeff * I) * (1 - 1 / (2**n_qpe)) / (2 * normalization_coeff)\n",
    "    )\n",
    "\n",
    "else:\n",
    "    Mnew = M"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f26c0aa2-2866-4f26-9d45-919e007c1b27",
   "metadata": {},
   "source": [
    "## 4. Building the Quantum Model\n",
    "\n",
    "Create a quantum model of the QPE algorithm using the Classiq platform with your desired constraints and preferences.\n",
    "\n",
    "There are generally two methods for inserting the matrix into the QFT: unitary implementation, which is exact but long; and exponentiation, which is approximated but shorter in depth. Choose the parameter `IS_EXACT` to indicate the chosen method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "67d1f58d-bb01-4041-bbf4-fdc52f30cf7f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:12.243903Z",
     "iopub.status.busy": "2024-05-07T14:54:12.243429Z",
     "iopub.status.idle": "2024-05-07T14:54:12.334591Z",
     "shell.execute_reply": "2024-05-07T14:54:12.333840Z"
    }
   },
   "outputs": [],
   "source": [
    "import scipy\n",
    "\n",
    "from classiq import Output, QNum, allocate, create_model, if_, prepare_amplitudes\n",
    "\n",
    "IS_EXACT = True\n",
    "\n",
    "my_amp = (\n",
    "    int_vec / np.linalg.norm(int_vec)\n",
    ").tolist()  # amplitude is given by the eignevector\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def main(phase_result: Output[QNum[n_qpe, False, n_qpe]]) -> None:\n",
    "    state = QArray(\"state\")\n",
    "    prepare_amplitudes(my_amp, 0.0, state)\n",
    "    phase_array = QArray(\"phase_array\")\n",
    "    my_qpe_flexible(\n",
    "        unitary_with_power=lambda pw, target: if_(\n",
    "            condition=IS_EXACT,\n",
    "            then=lambda: unitary_with_power_logic(\n",
    "                matrix=scipy.linalg.expm(1j * 2 * np.pi * Mnew).tolist(),\n",
    "                pw=pw,\n",
    "                target=target,\n",
    "            ),\n",
    "            else_=lambda: suzuki_trotter1_with_power_logic(\n",
    "                hamiltonian=pauli_operator_to_hamiltonian(pauli_ops),\n",
    "                pw=pw,\n",
    "                r0=2,\n",
    "                reps_scaling_factor=1.8,\n",
    "                evolution_coefficient=-2 * np.pi,\n",
    "                target=target,\n",
    "            ),\n",
    "        ),\n",
    "        precision=n_qpe,\n",
    "        state=state,\n",
    "        phase=phase_array,\n",
    "    )\n",
    "    bind(phase_array, phase_result)\n",
    "\n",
    "\n",
    "qmod = create_model(main)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bca846b9-b3c8-42d1-b3a8-5cb0c3f3b074",
   "metadata": {},
   "source": [
    "Set execution preferences:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "c8b701e3-1393-427f-b705-662fcfaf8728",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:12.340767Z",
     "iopub.status.busy": "2024-05-07T14:54:12.339097Z",
     "iopub.status.idle": "2024-05-07T14:54:12.381320Z",
     "shell.execute_reply": "2024-05-07T14:54:12.380565Z"
    }
   },
   "outputs": [],
   "source": [
    "from classiq import set_execution_preferences\n",
    "from classiq.execution import ExecutionPreferences\n",
    "\n",
    "num_shots = 10000\n",
    "qmod = set_execution_preferences(qmod, ExecutionPreferences(num_shots=num_shots))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "ca339f0c-4ffe-4130-8541-5621325e0c84",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:12.386254Z",
     "iopub.status.busy": "2024-05-07T14:54:12.385075Z",
     "iopub.status.idle": "2024-05-07T14:54:12.416936Z",
     "shell.execute_reply": "2024-05-07T14:54:12.416272Z"
    }
   },
   "outputs": [],
   "source": [
    "from classiq import write_qmod\n",
    "\n",
    "write_qmod(qmod, \"qpe_for_matrix\", decimal_precision=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f646cdf-c78f-4d28-822e-7eda986ad7cc",
   "metadata": {},
   "source": [
    "Synthesize the circuit and display it with the analyzer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "0ccdea1b-c255-42ab-83f2-9e36b561f94e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:12.421846Z",
     "iopub.status.busy": "2024-05-07T14:54:12.420689Z",
     "iopub.status.idle": "2024-05-07T14:54:19.781604Z",
     "shell.execute_reply": "2024-05-07T14:54:19.780782Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Opening: https://platform.classiq.io/circuit/11f3b909-2761-49cb-b3a4-1112d8036df6?version=0.41.0.dev39%2B79c8fd0855\n"
     ]
    }
   ],
   "source": [
    "from classiq import show, synthesize\n",
    "\n",
    "qprog = synthesize(qmod)\n",
    "show(qprog)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99e84bb8-90b8-4c40-9c04-51e18d7b86c9",
   "metadata": {},
   "source": [
    "# 5. Measuring and Analyzing the Generated Circuit \n",
    "\n",
    "Execute the circuit and analyze the results obtained from the quantum program, in comparison to the expected classical ones."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a9ad253-ac88-4c5e-a251-c40b2bd44bdb",
   "metadata": {},
   "source": [
    "### 5.1. Run the Circuit \n",
    "\n",
    "Send the circuit for execution by a chosen backend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "7a6b833b-ec74-4432-9bf5-0ce351149761",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:19.785625Z",
     "iopub.status.busy": "2024-05-07T14:54:19.784930Z",
     "iopub.status.idle": "2024-05-07T14:54:21.272353Z",
     "shell.execute_reply": "2024-05-07T14:54:21.271541Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from classiq import execute\n",
    "\n",
    "results = execute(qprog).result()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "5717e0dd-0354-4f88-be2c-88442a81ead0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:21.277738Z",
     "iopub.status.busy": "2024-05-07T14:54:21.276403Z",
     "iopub.status.idle": "2024-05-07T14:54:21.281650Z",
     "shell.execute_reply": "2024-05-07T14:54:21.280949Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from classiq.execution import ExecutionDetails\n",
    "\n",
    "results = results[0].value"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eddd2a37-9fe4-4646-8ab3-0c29f24ef5d1",
   "metadata": {},
   "source": [
    "Choose the number of eigenvalues to extract from the poll of results. The `number_of_solutions` value determines how many results from `qpe_results` are analyzed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "6ac801fd-5933-4b84-8a8b-6c520ba77842",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:21.286351Z",
     "iopub.status.busy": "2024-05-07T14:54:21.285175Z",
     "iopub.status.idle": "2024-05-07T14:54:21.289973Z",
     "shell.execute_reply": "2024-05-07T14:54:21.289304Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "number_of_solutions = 2  # number of phases sought"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c3b307b-a3b8-4103-a47f-f0fc3f1e8009",
   "metadata": {},
   "source": [
    "### 5.2. Translate into Eigenvalues (Phases)\n",
    "\n",
    "Here, the value in the `results` vector is translated from a binary number into a full solution for the eigenvalues."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ef03ea1-4ae0-4ce3-8b7c-aec30dcbba4d",
   "metadata": {},
   "source": [
    "Initially, use the parsed results to obtain the phases of the normalized matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "87718e4b-3ccb-470a-850b-0e3534396b28",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:21.294582Z",
     "iopub.status.busy": "2024-05-07T14:54:21.293405Z",
     "iopub.status.idle": "2024-05-07T14:54:21.312159Z",
     "shell.execute_reply": "2024-05-07T14:54:21.311469Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Your decimal solutions are\n",
      "[0.6875, 0.58203125]\n"
     ]
    }
   ],
   "source": [
    "dec_sol_vec = [\n",
    "    sampled_state.state[\"phase_result\"]\n",
    "    for sampled_state in results.parsed_counts[:number_of_solutions]\n",
    "]\n",
    "\n",
    "print(\"Your decimal solutions are\", dec_sol_vec, sep=\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cb35ab6-22d0-4ea3-8b6d-4bca00ca5c51",
   "metadata": {},
   "source": [
    "Then these decimal values are mapped back into the original values; i.e., renormalized into the original span."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "1409b9d6-4681-4eb9-9dc5-d9bca9c6ed1e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:21.316973Z",
     "iopub.status.busy": "2024-05-07T14:54:21.315786Z",
     "iopub.status.idle": "2024-05-07T14:54:21.321230Z",
     "shell.execute_reply": "2024-05-07T14:54:21.320573Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# renormalize into the \"real\" solution -\n",
    "if normalize:\n",
    "    solution = [\n",
    "        ((value * 2 * normalization_coeff / (1 - (1 / 2**n_qpe))) - normalization_coeff)\n",
    "        for value in dec_sol_vec\n",
    "    ]\n",
    "else:\n",
    "    solution = dec_sol_vec"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29a119e5-6883-43f8-abef-cb91ed934644",
   "metadata": {},
   "source": [
    "These are the results of the phases (matrix eigenvalues):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "2f8c746d-3659-4793-99e8-aec8250cecec",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:21.325798Z",
     "iopub.status.busy": "2024-05-07T14:54:21.324469Z",
     "iopub.status.idle": "2024-05-07T14:54:21.330800Z",
     "shell.execute_reply": "2024-05-07T14:54:21.330157Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.8935902927058823, 0.39612765552941154]\n"
     ]
    }
   ],
   "source": [
    "print(solution)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec60ea88-ef3a-41cc-bda7-945abe330ae3",
   "metadata": {},
   "source": [
    "These are the results, including the error contributed from the resolution (the number of qubits participating in the QPE):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "79f1c1ba-61fd-438e-a3d5-7606ea95e70f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:21.335165Z",
     "iopub.status.busy": "2024-05-07T14:54:21.334061Z",
     "iopub.status.idle": "2024-05-07T14:54:21.341661Z",
     "shell.execute_reply": "2024-05-07T14:54:21.340945Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the resolution of results is 0.018424542117647057\n",
      "the solutions are between 0.8751657505882352 and 0.9120148348235294\n",
      "the solutions are between 0.3777031134117645 and 0.4145521976470586\n"
     ]
    }
   ],
   "source": [
    "if normalize:\n",
    "    energy_resolution = (\n",
    "        (1 / (2**n_qpe)) * 2 * normalization_coeff / (1 - (1 / 2**n_qpe))\n",
    "    )\n",
    "else:\n",
    "    energy_resolution = 1 / (2**n_qpe)\n",
    "\n",
    "print(\"the resolution of results is\", energy_resolution)\n",
    "\n",
    "for sol in solution:\n",
    "    print(\n",
    "        \"the solutions are between\",\n",
    "        sol - energy_resolution,\n",
    "        \"and\",\n",
    "        sol + energy_resolution,\n",
    "    )\n",
    "\n",
    "    ### if zero or exceeds the normalization range, need to add conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51b1e9b4-d004-493e-8278-77c1b3cf0301",
   "metadata": {},
   "source": [
    "### 5.3. Compare to Exact Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "4032aaa8-d124-4c76-bea1-646e3162b02a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:21.345794Z",
     "iopub.status.busy": "2024-05-07T14:54:21.344679Z",
     "iopub.status.idle": "2024-05-07T14:54:21.352079Z",
     "shell.execute_reply": "2024-05-07T14:54:21.351463Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the eigenvalues are [0.9 0.4 0.1 0.2]\n",
      "the eigenvectors are\n",
      "[[ 0.51510515  0.41480695  0.5588446  -0.50029451]\n",
      " [ 0.61747259 -0.30596016 -0.64233734 -0.3354381 ]\n",
      " [ 0.58498122  0.11134965  0.09256217  0.79801659]\n",
      " [-0.10578874  0.8496616  -0.51626321  0.01887348]]\n"
     ]
    }
   ],
   "source": [
    "w, v = LA.eig(M)\n",
    "\n",
    "print(\"the eigenvalues are\", w)\n",
    "print(\"the eigenvectors are\", v, sep=\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb8793fc-6ef6-4ca7-9654-7d1267caa6e0",
   "metadata": {},
   "source": [
    "### 5.4. Find the Solution's Histogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "4b35f44e-a0a1-449b-9015-bf6af2ad53ea",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:21.356703Z",
     "iopub.status.busy": "2024-05-07T14:54:21.355583Z",
     "iopub.status.idle": "2024-05-07T14:54:21.364508Z",
     "shell.execute_reply": "2024-05-07T14:54:21.363833Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "energy_vec = []\n",
    "energy_prob = []\n",
    "\n",
    "for sampled_state in results.parsed_counts:\n",
    "    temp = sampled_state.state[\"phase_result\"]\n",
    "    if normalize:\n",
    "        temp2 = (\n",
    "            temp * 2 * normalization_coeff / (1 - (1 / 2**n_qpe))\n",
    "        ) - normalization_coeff\n",
    "    else:\n",
    "        temp2 = temp\n",
    "    energy_vec.append(temp2)\n",
    "    energy_prob.append(sampled_state.shots / num_shots)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "463a3503-dee1-4fb0-bdf7-e698aafd0b9f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:54:21.369146Z",
     "iopub.status.busy": "2024-05-07T14:54:21.368053Z",
     "iopub.status.idle": "2024-05-07T14:54:21.615666Z",
     "shell.execute_reply": "2024-05-07T14:54:21.614769Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAe5UlEQVR4nO3df2xV933w8c+1I2x+OvA4MYSYmSSrsgoCLRgXutCm8sqkrFNGVnk/NDyry6Y1QYq8SQM24XVVY7pGHVLgCWmeNZM6pfG6hyRq14dpc/NjWeiI4KFN2iWPkoaZAjZYTa+ZQbjyvc8fFCcG2/E113x97ddLukpycs693+N7fO/b555zbiafz+cDACCRstQDAABmNjECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJXZd6AOORy+Xi5MmTMX/+/MhkMqmHAwCMQz6fj7Nnz8ZNN90UZWWj7/8oiRg5efJk1NbWph4GADABx48fj5tvvnnU/18SMTJ//vyIuLgyCxYsSDwaAGA8+vr6ora2duh9fDQlESOXPppZsGCBGAGAEvN+h1g4gBUASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRgGniVPZ8vPxWb5zKnk89FChISXw3DQBj63ilK7bvfzVy+YiyTET75pXRVL8s9bBgXOwZAShxp7Lnh0IkIiKXj9ix/zV7SCgZYgSgxL3d2z8UIpcM5vNxrPdcmgFBgcQIQIlbXj03yi77hvbyTCbqquekGRAUSIwAlLglVbOjffPKKM9cLJLyTCYe2rwillTNTjwyGB8HsAJMA031y2LjB26IY73noq56jhChpIgRgGliSdVsEUJJ8jENAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFITipG9e/dGXV1dVFZWRkNDQxw6dGhcyz311FORyWTinnvumcjDAgDTUMEx0tHREa2trdHW1hZHjhyJVatWxaZNm+L06dNjLnfs2LH40z/907jzzjsnPFgAYPopOEa+/OUvx3333RctLS3xwQ9+MPbt2xdz5syJr371q6MuMzg4GL/7u78bn/vc5+KWW265qgEDANNLQTEyMDAQhw8fjsbGxnfvoKwsGhsb4+DBg6Mu91d/9Vdx4403xmc+85lxPc6FCxeir69v2A0AmJ4KipHe3t4YHByMmpqaYdNramqiu7t7xGVeeuml+Nu//dt4/PHHx/047e3tUVVVNXSrra0tZJgAQAmZ1LNpzp49G7/3e78Xjz/+eFRXV497ue3bt0c2mx26HT9+fBJHCQCkdF0hM1dXV0d5eXn09PQMm97T0xOLFy++Yv633norjh07Fp/61KeGpuVyuYsPfN118cYbb8Stt956xXIVFRVRUVFRyNAAgBJV0J6RWbNmxZo1a6Kzs3NoWi6Xi87Ozli/fv0V899+++3x6quvxtGjR4duv/7rvx533XVXHD161McvAEBhe0YiIlpbW6O5uTnWrl0b69ati927d0d/f3+0tLRERMSWLVti6dKl0d7eHpWVlbFixYphy19//fUREVdMBwBmpoJjpKmpKc6cORM7d+6M7u7uWL16dRw4cGDooNaurq4oK3NhVwBgfDL5fD6fehDvp6+vL6qqqiKbzcaCBQtSDwcAGIfxvn/bhQEAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKQmFCN79+6Nurq6qKysjIaGhjh06NCo8+7fvz/Wrl0b119/fcydOzdWr14dX/va1yY8YABgeik4Rjo6OqK1tTXa2triyJEjsWrVqti0aVOcPn16xPkXLVoUf/7nfx4HDx6M73//+9HS0hItLS3xz//8z1c9eACg9GXy+Xy+kAUaGhqivr4+9uzZExERuVwuamtrY+vWrbFt27Zx3ceHP/zhuPvuu+Pzn//8uObv6+uLqqqqyGazsWDBgkKGCwAkMt7374L2jAwMDMThw4ejsbHx3TsoK4vGxsY4ePDg+y6fz+ejs7Mz3njjjdi4ceOo8124cCH6+vqG3QCA6amgGOnt7Y3BwcGoqakZNr2mpia6u7tHXS6bzca8efNi1qxZcffdd8cjjzwSv/IrvzLq/O3t7VFVVTV0q62tLWSYAEAJuSZn08yfPz+OHj0ar7zySnzhC1+I1tbWeP7550edf/v27ZHNZodux48fvxbDBAASuK6Qmaurq6O8vDx6enqGTe/p6YnFixePulxZWVncdtttERGxevXq+M///M9ob2+Pj3/84yPOX1FRERUVFYUMDQAoUQXtGZk1a1asWbMmOjs7h6blcrno7OyM9evXj/t+crlcXLhwoZCHBgCmqYL2jEREtLa2RnNzc6xduzbWrVsXu3fvjv7+/mhpaYmIiC1btsTSpUujvb09Ii4e/7F27dq49dZb48KFC/Htb387vva1r8Wjjz5a3DUBAEpSwTHS1NQUZ86ciZ07d0Z3d3esXr06Dhw4MHRQa1dXV5SVvbvDpb+/Pz772c/Gj3/845g9e3bcfvvt8fd///fR1NRUvLUAKIJT2fPxdm9/LK+eG0uqZqceDswYBV9nJAXXGQEmW8crXbF9/6uRy0eUZSLaN6+MpvplqYcFJW1SrjMCMB2dyp4fCpGIiFw+Ysf+1+JU9nzagcEMIUaAGe/t3v6hELlkMJ+PY73n0gwIZhgxAsx4y6vnRllm+LTyTCbqquekGRDMMGIEmPGWVM2O9s0rozxzsUjKM5l4aPMKB7HCNVLw2TQA01FT/bLY+IEb4ljvuairniNE4BoSIwA/t6RqtgiBBHxMAwAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSmlCM7N27N+rq6qKysjIaGhri0KFDo877+OOPx5133hkLFy6MhQsXRmNj45jzAwAzS8Ex0tHREa2trdHW1hZHjhyJVatWxaZNm+L06dMjzv/888/Hb//2b8dzzz0XBw8ejNra2vjkJz8ZJ06cuOrBAwClL5PP5/OFLNDQ0BD19fWxZ8+eiIjI5XJRW1sbW7dujW3btr3v8oODg7Fw4cLYs2dPbNmyZVyP2dfXF1VVVZHNZmPBggWFDBcASGS8798F7RkZGBiIw4cPR2Nj47t3UFYWjY2NcfDgwXHdx7lz5+JnP/tZLFq0aNR5Lly4EH19fcNuAMD0VFCM9Pb2xuDgYNTU1AybXlNTE93d3eO6jz/7sz+Lm266aVjQXK69vT2qqqqGbrW1tYUMEwAoIdf0bJpdu3bFU089FU8//XRUVlaOOt/27dsjm80O3Y4fP34NRwkAXEvXFTJzdXV1lJeXR09Pz7DpPT09sXjx4jGXffjhh2PXrl3xr//6r3HHHXeMOW9FRUVUVFQUMjQAoEQVtGdk1qxZsWbNmujs7ByalsvlorOzM9avXz/qcn/9138dn//85+PAgQOxdu3aiY8WAJh2CtozEhHR2toazc3NsXbt2li3bl3s3r07+vv7o6WlJSIitmzZEkuXLo329vaIiPjiF78YO3fujCeffDLq6uqGji2ZN29ezJs3r4irAgCUooJjpKmpKc6cORM7d+6M7u7uWL16dRw4cGDooNaurq4oK3t3h8ujjz4aAwMD8Zu/+ZvD7qetrS3+8i//8upGDwCUvIKvM5KC64wAQOmZlOuMAAAUmxgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEYJo5lT0fL7/VG6ey51MPBcblutQDAKB4Ol7piu37X41cPqIsE9G+eWU01S9LPSwYkz0jANPEqez5oRCJiMjlI3bsf80eEqY8MQIwTbzd2z8UIpcM5vNxrPdcmgHBOIkRgGliefXcKMsMn1aeyURd9Zw0A4JxEiMA08SSqtnRvnlllGcuFkl5JhMPbV4RS6pmJx4ZjM0BrADTSFP9stj4gRviWO+5qKueI0QoCWIEYJpZUjVbhFBSfEwDACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJTShG9u7dG3V1dVFZWRkNDQ1x6NChUef9wQ9+EPfee2/U1dVFJpOJ3bt3T3SsAMA0VHCMdHR0RGtra7S1tcWRI0di1apVsWnTpjh9+vSI8587dy5uueWW2LVrVyxevPiqBwwATC8Fx8iXv/zluO+++6KlpSU++MEPxr59+2LOnDnx1a9+dcT56+vr40tf+lL81m/9VlRUVFz1gAGA6aWgGBkYGIjDhw9HY2Pju3dQVhaNjY1x8ODBog3qwoUL0dfXN+wGAExPBcVIb29vDA4ORk1NzbDpNTU10d3dXbRBtbe3R1VV1dCttra2aPcNAEwtU/Jsmu3bt0c2mx26HT9+PPWQAIBJUtC39lZXV0d5eXn09PQMm97T01PUg1MrKiocXwIAM0RBe0ZmzZoVa9asic7OzqFpuVwuOjs7Y/369UUfHAAw/RW0ZyQiorW1NZqbm2Pt2rWxbt262L17d/T390dLS0tERGzZsiWWLl0a7e3tEXHxoNcf/vCHQ/9+4sSJOHr0aMybNy9uu+22Iq4KAFCKCo6RpqamOHPmTOzcuTO6u7tj9erVceDAgaGDWru6uqKs7N0dLidPnowPfehDQ//98MMPx8MPPxwf+9jH4vnnn7/6NQAASlomn8/nUw/i/fT19UVVVVVks9lYsGBB6uEAAOMw3vfvKXk2DQAwc4gRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBKHGnsufj5bd641T2fOqhwIRcl3oAAExcxytdsX3/q5HLR5RlIto3r4ym+mWphwUFsWcEoESdyp4fCpGIiFw+Ysf+1+whoeSIEYAS9XZv/1CIXDKYz8ex3nNpBgQTJEYAStTy6rlRlhk+rTyTibrqOWkGBBMkRgBK1JKq2dG+eWWUZy4WSXkmEw9tXhFLqmYnHhkUxgGsACWsqX5ZbPzADXGs91zUVc8RIpQkMQJQ4pZUzRYhlDQf0wBcxnU74NqyZwTgPabTdTtOZc/H2739sbx6rj0nTGliBODnRrtux8YP3FByb+bTKaqY/nxMA/Bz0+W6HS6GRqkRIwA/N12u2zFdooqZQ4wAU1KKg0iny3U7pktUMXM4ZgSYclIe7zAdrttxKap27H8tBvP5ko0qZo5MPp/Pv/9safX19UVVVVVks9lYsGBB6uEAk+hU9nx8dNd3hn3MUJ7JxEvb7vJmOoKxzpg5lT1f0lFF6Rvv+7c9I8CUMtbxDpP1hlqqp8C+3x4kF0OjVIgRYEq5dLzD5XtGJut4h1I9BbaQ05BLNbaYORzACkwp1/Ig0lI+BXa8Z8x0vNIVH931nfidx/8jPrrrO9HxStc1HCWMjz0jwJRzrQ4iTfGRULHMnVUemYh47/Av34M0nS7ixvQmRoAp6Voc7zDWR0JT9aONU9nz8cRLb8fj//b2FSFy+R6kUo4tZpYJfUyzd+/eqKuri8rKymhoaIhDhw6NOf83vvGNuP3226OysjJWrlwZ3/72tyc0WIBi+4NfXj70QnjpDf3F/3dm6KONDe3fiYf+6YfJP7o5lT0fD/3TD2ND+3fiK5eFSFkmYv9n119xrMtI1xuJiPjm90/E946/48sAmTIKPrW3o6MjtmzZEvv27YuGhobYvXt3fOMb34g33ngjbrzxxivmf/nll2Pjxo3R3t4ev/ZrvxZPPvlkfPGLX4wjR47EihUrxvWYk3VqbyF/+Uzkr6Sr/ctqpOWL/ddasdZrsn+WxXjcyRjTZCnG81KM7e/wf70T+Xw+1tYtGvU+J/Lcz51VHv0Dg1f8c3n13IiIEddjpHlGup+un5yLTCYTtQtnj/gYc2eVx/F3zsfLb/bG1w8dj3xcfDP/g1++JVp+uS4i4opTiyMiMhFx353L4+47lhT0WGPNM5Gxjubr930k1t/6P66Y/tgLb0X7/3l91OXeu17H3zkf+Xw+li2aM6H1ee+Yx7qfsbafsbaxkZYfj9G20fc7LXq0xxzPcqlfS67Fa/J4jPf9u+AYaWhoiPr6+tizZ09ERORyuaitrY2tW7fGtm3brpi/qakp+vv741vf+tbQtI985COxevXq2LdvX1FXphCFHEE/kaPtr/YI/ZGWj4iiHvVfrPUqZFwT/blM5s9jqp1NUYzn5Tc+tDSe/r8nrmr72/a/Xx1688tExOYPX3mfERN77kdz6Y/4S4Hw3vUYaZ5iunQtk7d7++N3Hv+PIt/75BrrOiwvv9U75dZntO1npO32vfNcvvxEXzua6peN+Xs20rY6nm1+qryWTPb7WyEmJUYGBgZizpw58Y//+I9xzz33DE1vbm6On/70p/Hss89escyyZcuitbU1HnzwwaFpbW1t8cwzz8T3vve9oq7MeBVyUaWJXIDpai/aNNLyZZmIfP7Kg9UmeiGoYq1XWUTECJ+5F+tnWYzHHctUu8BWsZ6Xy13t9jeSQrbJ8d5nal+/7yNRVz2nJMZ6SVlEtN87+hvIVP3Zj7T9FDLPRF87yjOZ2P/Z9fEb//PlEX/PIkbeMzbaeMZaLsVryWS/vxVqvO/fBR0z0tvbG4ODg1FTUzNsek1NTXR3d4+4THd3d0HzR0RcuHAh+vr6ht2KqZAvkZrIF05d7ZdUjbR8boRfyKv54qtirVcuYlJ/lsV43ELvO+UXihXrebnc1W5/IylkmxzvfaZ06cDVS6cWT/XrHmQi4g/vvCX+ffsnxvxLdqquz0jbTyHzTPS1YzCfj1eOvTPq79lY2+pY2/xUeS2Z7Pe3yTIlz6Zpb2+Pz33uc5N2/4VcVGkiF2C62os2jbT8aEU+0QtBFWu9RttDUayfZTEet9D7TvmFYsV6Xi53tdvfSArZJsd7n6mURQw7E+XSqcVPvHQsHv+3HxX9I6GrcfEYj4vHt4z3r9f3rs//eulHU+J5KMaekYm8dpRnMlFft3DM37PRttX32+anwmvJZL+/TZaCYrm6ujrKy8ujp6dn2PSenp5YvHjxiMssXry4oPkjIrZv3x7ZbHbodvz48UKG+b4KuajSRC7AdLUXbRpp+fbNK2PXvcW7EFSx1qv93pWT+rMsxuMW++cwmYr1vNz74aVXvf1l3nMWRibiivssZJu8fIyjyUQMPe7l6zE0T+bd40au1lh7F5ZUzY4dd/9SvLz9E/GHd94y4lkp19Klsb68/ROx4+5fmtC2vuPuX4p/3/aJ+Pp9H4ln79+QbL1G237ebxt77/ITfe14aPOKWFW7cNTfs9G21ffb5qfKa8lkv79NlgkdwLpu3bp45JFHIuLiAazLli2LBx54YNQDWM+dOxff/OY3h6Zt2LAh7rjjjqQHsEYU9iVSE/nCqav9kqqRli/2F18Va70m+2dZjMedjDFNlmI8L8XY/o781zuRz0esqVs46n1O5LmfM6sszg3krvjnpb/IRlqPkeYZ6X6O/+R8ZDIRNy+cPeJjzJlVFj9+5/yw9Srk51voY401z2SNtRDvXa9Lj1W7aGLrc/mYR7ufsbafsbaxkZYvZB0vX26sbXesxxzPcqlfS67Fa/J4TNrZNB0dHdHc3ByPPfZYrFu3Lnbv3h3/8A//EK+//nrU1NTEli1bYunSpdHe3h4RF0/t/djHPha7du2Ku+++O5566ql46KGHpsSpvQDA5Jm0b+1tamqKM2fOxM6dO6O7uztWr14dBw4cGDpItaurK8rK3v30Z8OGDfHkk0/GX/zFX8SOHTviF3/xF+OZZ54Zd4gAANNbwXtGUrBnBABKz6Sc2gsAUGxiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSBV8OPoVLF4nt6+tLPBIAYLwuvW+/38XeSyJGzp49GxERtbW1iUcCABTq7NmzUVVVNer/L4nvpsnlcnHy5MmYP39+ZDKZ1MOZkL6+vqitrY3jx4/7fp2EPA/peQ6mBs9DejPhOcjn83H27Nm46aabhn2J7uVKYs9IWVlZ3HzzzamHURQLFiyYthtdKfE8pOc5mBo8D+lN9+dgrD0ilziAFQBISowAAEmJkWukoqIi2traoqKiIvVQZjTPQ3qeg6nB85Ce5+BdJXEAKwAwfdkzAgAkJUYAgKTECACQlBgBAJISI9fYsWPH4jOf+UwsX748Zs+eHbfeemu0tbXFwMBA6qHNOF/4whdiw4YNMWfOnLj++utTD2fG2Lt3b9TV1UVlZWU0NDTEoUOHUg9pRnnxxRfjU5/6VNx0002RyWTimWeeST2kGae9vT3q6+tj/vz5ceONN8Y999wTb7zxRuphJSVGrrHXX389crlcPPbYY/GDH/wg/uZv/ib27dsXO3bsSD20GWdgYCA+/elPxx//8R+nHsqM0dHREa2trdHW1hZHjhyJVatWxaZNm+L06dOphzZj9Pf3x6pVq2Lv3r2phzJjvfDCC3H//ffHd7/73fiXf/mX+NnPfhaf/OQno7+/P/XQknFq7xTwpS99KR599NH40Y9+lHooM9Lf/d3fxYMPPhg//elPUw9l2mtoaIj6+vrYs2dPRFz83qna2trYunVrbNu2LfHoZp5MJhNPP/103HPPPamHMqOdOXMmbrzxxnjhhRdi48aNqYeThD0jU0A2m41FixalHgZMqoGBgTh8+HA0NjYOTSsrK4vGxsY4ePBgwpFBWtlsNiJiRr8PiJHE3nzzzXjkkUfij/7oj1IPBSZVb29vDA4ORk1NzbDpNTU10d3dnWhUkFYul4sHH3wwPvrRj8aKFStSDycZMVIk27Zti0wmM+bt9ddfH7bMiRMn4ld/9Vfj05/+dNx3332JRj69TOR5AEjl/vvvj9deey2eeuqp1ENJ6rrUA5gu/uRP/iR+//d/f8x5brnllqF/P3nyZNx1112xYcOG+MpXvjLJo5s5Cn0euHaqq6ujvLw8enp6hk3v6emJxYsXJxoVpPPAAw/Et771rXjxxRfj5ptvTj2cpMRIkdxwww1xww03jGveEydOxF133RVr1qyJJ554IsrK7KAqlkKeB66tWbNmxZo1a6Kzs3PogMlcLhednZ3xwAMPpB0cXEP5fD62bt0aTz/9dDz//POxfPny1ENKToxcYydOnIiPf/zj8Qu/8Avx8MMPx5kzZ4b+n78Or62urq74yU9+El1dXTE4OBhHjx6NiIjbbrst5s2bl3Zw01Rra2s0NzfH2rVrY926dbF79+7o7++PlpaW1EObMf77v/873nzzzaH/fvvtt+Po0aOxaNGiWLZsWcKRzRz3339/PPnkk/Hss8/G/Pnzh46ZqqqqitmzZyceXSJ5rqknnngiHxEj3ri2mpubR3wennvuudRDm9YeeeSR/LJly/KzZs3Kr1u3Lv/d73439ZBmlOeee27E7b65uTn10GaM0d4DnnjiidRDS8Z1RgCApBysAAAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCS+v+EJowS3g2+4wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(energy_vec, energy_prob, \".\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f2e9b12-07fd-4105-beee-620ba6a00c20",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "<a name='NC'>[1]</a>: [Michael A. Nielsen and Isaac L. Chuang. 2011. Quantum Computation and Quantum Information: 10th Anniversary Edition, Cambridge University Press, New York, NY, USA.\n",
    "](http://mmrc.amss.cas.cn/tlb/201702/W020170224608149940643.pdf)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
